##### Replication for Weidong, Zhang. "The importance of money and connections: 
# Explaining international status from UNGA draft sponsorship networks" (1)
#### Feb, 2025


### last updated  2/13, 2025


setwd("~/Desktop/submission/FDI Networks/International Interaction/Third Revision-1:31:25/Replication Materials")
# Note: change it to your own working directory
# Estimated simulation time: 3 hours


# The analysis is conducted by using R version 4.1.1 (2021-08-10) -- "Kick Things"


library("writexl")
library(plyr)
library(dplyr)
library(tidyverse)
library(haven)
library(countrycode)
library(igraph)
library(ergm)
library(btergm) 
library(readxl)
library(tibble)
library(ggplot2)



# Descriptive statistics Table 1 ====
sp2009 <- read.csv("Edge2009.csv",row.names=1)

sp2010 <- read.csv("Edge2010.csv",row.names=1)

sp2011 <- read.csv("Edge2011.csv",row.names=1)

sp2012 <- read.csv("Edge2012.csv",row.names=1)
sp2013 <- read.csv("Edge2013.csv",row.names=1)
sp2014 <- read.csv("Edge2014.csv",row.names=1)

sp2015 <- read.csv("Edge2015.csv",row.names=1)
sp2016 <- read.csv("Edge2016.csv",row.names=1)
sp2017 <- read.csv("Edge2017.csv",row.names=1)
sp2018 <- read.csv("Edge2018.csv",row.names=1)




g2009 <- graph.adjacency(as.matrix(sp2009), mode="directed",weighted=NULL)

g2010 <- graph.adjacency(as.matrix(sp2010), mode="directed",weighted=NULL)

g2011 <- graph.adjacency(as.matrix(sp2011), mode="directed",weighted=NULL)

g2012 <- graph.adjacency(as.matrix(sp2012), mode="directed",weighted=NULL)
g2013 <- graph.adjacency(as.matrix(sp2013), mode="directed",weighted=NULL)
g2014 <- graph.adjacency(as.matrix(sp2014), mode="directed",weighted=NULL)

g2015 <- graph.adjacency(as.matrix(sp2015), mode="directed",weighted=NULL)
g2016 <- graph.adjacency(as.matrix(sp2016), mode="directed",weighted=NULL)
g2017 <- graph.adjacency(as.matrix(sp2017), mode="directed",weighted=NULL)
g2018 <- graph.adjacency(as.matrix(sp2018), mode="directed",weighted=NULL)



# number of isolates
number_isolates <- c(sum(degree(g2009)==0), sum(degree(g2010)==0), sum(degree(g2011)==0), 
                     sum(degree(g2012)==0), sum(degree(g2013)==0), sum(degree(g2014)==0), 
                     sum(degree(g2015)==0), sum(degree(g2016)==0), sum(degree(g2017)==0),sum(degree(g2018)==0))
number_isolates



# number of ties
number_ties <- c(ecount(g2009), ecount(g2010), ecount(g2011), ecount(g2012), ecount(g2013), ecount(g2014), 
                 ecount(g2015), ecount(g2016), ecount(g2017),ecount(g2018))
number_ties

# density
density <- c(graph.density(g2009), graph.density(g2010), graph.density(g2011), graph.density(g2012), 
             graph.density(g2013), graph.density(g2014), graph.density(g2015), graph.density(g2016),
             graph.density(g2017),graph.density(g2018))
density

# proportion of symmetric dyads
sym_dyads <- c(1-(dyad.census(g2009)$asym/(vcount(g2009)*(vcount(g2009)-1)/2)),
               1-(dyad.census(g2010)$asym/(vcount(g2010)*(vcount(g2010)-1)/2)),
               1-(dyad.census(g2011)$asym/(vcount(g2011)*(vcount(g2011)-1)/2)),
               1-(dyad.census(g2012)$asym/(vcount(g2012)*(vcount(g2012)-1)/2)),
               1-(dyad.census(g2013)$asym/(vcount(g2013)*(vcount(g2013)-1)/2)),
               1-(dyad.census(g2014)$asym/(vcount(g2014)*(vcount(g2014)-1)/2)),
               1-(dyad.census(g2015)$asym/(vcount(g2015)*(vcount(g2015)-1)/2)),
               1-(dyad.census(g2016)$asym/(vcount(g2016)*(vcount(g2016)-1)/2)),
               1-(dyad.census(g2017)$asym/(vcount(g2017)*(vcount(g2017)-1)/2)),
               1-(dyad.census(g2018)$asym/(vcount(g2018)*(vcount(g2018)-1)/2)))
sym_dyads

# global clustering coefficient
clustering <- c(transitivity(g2009), transitivity(g2010), transitivity(g2011), transitivity(g2012), 
                transitivity(g2013), transitivity(g2014), transitivity(g2015), transitivity(g2016), 
                transitivity(g2017),transitivity(g2018)) 
clustering


# Prepare data for analysis
# For network in each year====
# Network 2018====
Edge2018 <- read.csv("Edge2018.csv", header=T, row.names=1, stringsAsFactors = FALSE)

# For node attributes (GDP, democracy, military expenditure) of 2018
load("vcov3.RData")

attri2018 <- subset(vcov3, year==2018,
                    select=state: log_military)


ally2018 <- read.csv("ally2018.csv", header=T, row.names=1)
ally2018[,1:196] <- lapply(ally2018[,1:196], as.numeric)

FDI2018 <- read.csv("FDI2018.csv", header=T, row.names=1)
FDI2018[,1:196] <- lapply(FDI2018[,1:196], as.numeric)

trade2018 <- read.csv("trade2018.csv", header=T, row.names=1)
trade2018[,1:196] <- lapply(trade2018[,1:196], as.numeric)


dis2018 <- read.csv("dis2018.csv", header=T, row.names=1)
dis2018[,1:196] <- lapply(dis2018[,1:196], as.numeric)

#set attributes

net2018 <- network(as.matrix(Edge2018), directed = T)

set.vertex.attribute(net2018,names(attri2018),attri2018)

set.network.attribute(net2018,"ally2018",as.matrix(ally2018))

set.network.attribute(net2018,"FDI2018",as.matrix(FDI2018))

set.network.attribute(net2018,"trade2018",as.matrix(trade2018))

set.network.attribute(net2018,"dis2018",as.matrix(dis2018))

# Network 2017====
Edge2017 <- read.csv("Edge2017.csv", header=T, row.names=1)

ally2017 <- read.csv("ally2017.csv", header=T, row.names=1)
ally2017[,1:196] <- lapply(ally2017[,1:196], as.numeric)

FDI2017 <- read.csv("FDI2017.csv", header=T, row.names=1)
FDI2017[,1:196] <- lapply(FDI2017[,1:196], as.numeric)

trade2017 <- read.csv("trade2017.csv", header=T, row.names=1)
trade2017[,1:196] <- lapply(trade2017[,1:196], as.numeric)


dis2017 <- read.csv("dis2017.csv", header=T, row.names=1)
dis2017[,1:196] <- lapply(dis2017[,1:196], as.numeric)


# For node attributes (GDP, democracy, military expenditure) of 2017
load("vcov3.RData")
attri2017 <- subset(vcov3, year==2017,
                    select=state: log_military)



net2017 <- network(as.matrix(Edge2017), directed = T)

set.vertex.attribute(net2017,names(attri2017),attri2017)

set.network.attribute(net2017,"ally2017",as.matrix(ally2017))

set.network.attribute(net2017,"FDI2017",as.matrix(FDI2017))

set.network.attribute(net2017,"trade2017",as.matrix(trade2017))

set.network.attribute(net2017,"dis2017",as.matrix(dis2017))


# Network 2016====
Edge2016 <- read.csv("Edge2016.csv", header=T, row.names=1)
ally2016 <- read.csv("ally2016.csv", header=T, row.names=1)
ally2016[,1:196] <- lapply(ally2016[,1:196], as.numeric)

FDI2016 <- read.csv("FDI2016.csv", header=T, row.names=1)
FDI2016[,1:196] <- lapply(FDI2016[,1:196], as.numeric)

trade2016 <- read.csv("trade2016.csv", header=T, row.names=1)
trade2016[,1:196] <- lapply(trade2016[,1:196], as.numeric)

dis2016 <- read.csv("dis2016.csv", header=T, row.names=1)
dis2016[,1:196] <- lapply(dis2016[,1:196], as.numeric)


# For node attributes (GDP, democracy, military expenditure) of 2016
load("vcov3.RData")
attri2016 <- subset(vcov3, year==2016,
                    select=state: log_military)

net2016 <- network(as.matrix(Edge2016), directed = T)

set.vertex.attribute(net2016,names(attri2016),attri2016)

set.network.attribute(net2016,"ally2016",as.matrix(ally2016))

set.network.attribute(net2016,"FDI2016",as.matrix(FDI2016))

set.network.attribute(net2016,"trade2016",as.matrix(trade2016))

set.network.attribute(net2016,"dis2016",as.matrix(dis2016))

# Network 2015====
Edge2015 <- read.csv("Edge2015.csv", header=T, row.names=1)
ally2015 <- read.csv("ally2015.csv", header=T, row.names=1)
ally2015[,1:196] <- lapply(ally2015[,1:196], as.numeric)

FDI2015 <- read.csv("FDI2015.csv", header=T, row.names=1)
FDI2015[,1:196] <- lapply(FDI2015[,1:196], as.numeric)

trade2015 <- read.csv("trade2015.csv", header=T, row.names=1)
trade2015[,1:196] <- lapply(trade2015[,1:196], as.numeric)

dis2015 <- read.csv("dis2015.csv", header=T, row.names=1)
dis2015[,1:196] <- lapply(dis2015[,1:196], as.numeric)


# For node attributes (GDP, democracy, military expenditure) of 2015
load("vcov3.RData")
attri2015 <- subset(vcov3, year==2015,
                    select=state: log_military)

net2015 <- network(as.matrix(Edge2015), directed = T)

set.vertex.attribute(net2015,names(attri2015),attri2015)

set.network.attribute(net2015,"ally2015",as.matrix(ally2015))

set.network.attribute(net2015,"FDI2015",as.matrix(FDI2015))

set.network.attribute(net2015,"trade2015",as.matrix(trade2015))

set.network.attribute(net2015,"dis2015",as.matrix(dis2015))


# Network 2014====
Edge2014 <- read.csv("Edge2014.csv", header=T, row.names=1)
ally2014 <- read.csv("ally2014.csv", header=T, row.names=1)
ally2014[,1:196] <- lapply(ally2014[,1:196], as.numeric)

FDI2014 <- read.csv("FDI2014.csv", header=T, row.names=1)
FDI2014[,1:196] <- lapply(FDI2014[,1:196], as.numeric)

trade2014 <- read.csv("trade2014.csv", header=T, row.names=1)
trade2014[,1:196] <- lapply(trade2014[,1:196], as.numeric)

dis2014 <- read.csv("dis2014.csv", header=T, row.names=1)
dis2014[,1:196] <- lapply(dis2014[,1:196], as.numeric)


# For node attributes (GDP, democracy, military expenditure) of 2014
load("vcov3.RData")
attri2014 <- subset(vcov3, year==2014,
                    select=state: log_military)

net2014 <- network(as.matrix(Edge2014), directed = T)

set.vertex.attribute(net2014,names(attri2014),attri2014)

set.network.attribute(net2014,"ally2014",as.matrix(ally2014))

set.network.attribute(net2014,"FDI2014",as.matrix(FDI2014))

set.network.attribute(net2014,"trade2014",as.matrix(trade2014))

set.network.attribute(net2014,"dis2014",as.matrix(dis2014))

# Network 2013====
Edge2013 <- read.csv("Edge2013.csv", header=T, row.names=1)
ally2013 <- read.csv("ally2013.csv", header=T, row.names=1)
ally2013[,1:196] <- lapply(ally2013[,1:196], as.numeric)

FDI2013 <- read.csv("FDI2013.csv", header=T, row.names=1)
FDI2013[,1:196] <- lapply(FDI2013[,1:196], as.numeric)
trade2013 <- read.csv("trade2013.csv", header=T, row.names=1)
trade2013[,1:196] <- lapply(trade2013[,1:196], as.numeric)

dis2013 <- read.csv("dis2013.csv", header=T, row.names=1)
dis2013[,1:196] <- lapply(dis2013[,1:196], as.numeric)


# For node attributes (GDP, democracy, military expenditure) of 2013
load("vcov3.RData")
attri2013 <- subset(vcov3, year==2013,
                    select=state: log_military)

net2013 <- network(as.matrix(Edge2013), directed = T)

set.vertex.attribute(net2013,names(attri2013),attri2013)

set.network.attribute(net2013,"ally2013",as.matrix(ally2013))

set.network.attribute(net2013,"FDI2013",as.matrix(FDI2013))
set.network.attribute(net2013,"trade2013",as.matrix(trade2013))

set.network.attribute(net2013,"dis2013",as.matrix(dis2013))


# Network 2012====
Edge2012 <- read.csv("Edge2012.csv", header=T, row.names=1)
ally2012 <- read.csv("ally2012.csv", header=T, row.names=1)
ally2012[,1:196] <- lapply(ally2012[,1:196], as.numeric)

FDI2012 <- read.csv("FDI2012.csv", header=T, row.names=1)
FDI2012[,1:196] <- lapply(FDI2012[,1:196], as.numeric)
trade2012 <- read.csv("trade2012.csv", header=T, row.names=1)
trade2012[,1:196] <- lapply(trade2012[,1:196], as.numeric)

dis2012 <- read.csv("dis2012.csv", header=T, row.names=1)
dis2012[,1:196] <- lapply(dis2012[,1:196], as.numeric)



# For node attributes (GDP, democracy, military expenditure) of 2012
load("vcov3.RData")
attri2012 <- subset(vcov3, year==2012,
                    select=state: log_military)

net2012 <- network(as.matrix(Edge2012), directed = T)

set.vertex.attribute(net2012,names(attri2012),attri2012)

set.network.attribute(net2012,"ally2012",as.matrix(ally2012))

set.network.attribute(net2012,"FDI2012",as.matrix(FDI2012))

set.network.attribute(net2012,"trade2012",as.matrix(trade2012))

set.network.attribute(net2012,"dis2012",as.matrix(dis2012))


# Network 2011====
Edge2011 <- read.csv("Edge2011.csv", header=T, row.names=1)
ally2011 <- read.csv("ally2011.csv", header=T, row.names=1)
ally2011[,1:196] <- lapply(ally2011[,1:196], as.numeric)

FDI2011 <- read.csv("FDI2011.csv", header=T, row.names=1)
FDI2011[,1:196] <- lapply(FDI2011[,1:196], as.numeric)
trade2011 <- read.csv("trade2011.csv", header=T, row.names=1)
trade2011[,1:196] <- lapply(trade2011[,1:196], as.numeric)

dis2011 <- read.csv("dis2011.csv", header=T, row.names=1)
dis2011[,1:196] <- lapply(dis2011[,1:196], as.numeric)


# For node attributes (GDP, democracy, military expenditure) of 2011
load("vcov3.RData")
attri2011 <- subset(vcov3, year==2011,
                    select=state: log_military)

net2011 <- network(as.matrix(Edge2011), directed = T)

set.vertex.attribute(net2011,names(attri2011),attri2011)

set.network.attribute(net2011,"ally2011",as.matrix(ally2011))

set.network.attribute(net2011,"FDI2011",as.matrix(FDI2011))
set.network.attribute(net2011,"trade2011",as.matrix(trade2011))

set.network.attribute(net2011,"dis2011",as.matrix(dis2011))


# Network 2010====
Edge2010 <- read.csv("Edge2010.csv", header=T, row.names=1)
ally2010 <- read.csv("ally2010.csv", header=T, row.names=1)
ally2010[,1:196] <- lapply(ally2010[,1:196], as.numeric)

FDI2010 <- read.csv("FDI2010.csv", header=T, row.names=1)
FDI2010[,1:196] <- lapply(FDI2010[,1:196], as.numeric)
trade2010 <- read.csv("trade2010.csv", header=T, row.names=1)
trade2010[,1:196] <- lapply(trade2010[,1:196], as.numeric)

dis2010 <- read.csv("dis2010.csv", header=T, row.names=1)
dis2010[,1:196] <- lapply(dis2010[,1:196], as.numeric)


# For node attributes (GDP, democracy, military expenditure) of 2010
load("vcov3.RData")
attri2010 <- subset(vcov3, year==2010,
                    select=state: log_military)

net2010 <- network(as.matrix(Edge2010), directed = T)

set.vertex.attribute(net2010,names(attri2010),attri2010)

set.network.attribute(net2010,"ally2010",as.matrix(ally2010))

set.network.attribute(net2010,"FDI2010",as.matrix(FDI2010))
set.network.attribute(net2010,"trade2010",as.matrix(trade2010))

set.network.attribute(net2010,"dis2010",as.matrix(dis2010))


# Network 2009====
Edge2009 <- read.csv("Edge2009.csv", header=T, row.names=1)
ally2009 <- read.csv("ally2009.csv", header=T, row.names=1)
ally2009[,1:196] <- lapply(ally2009[,1:196], as.numeric)

FDI2009 <- read.csv("FDI2009.csv", header=T, row.names=1)
FDI2009[,1:196] <- lapply(FDI2009[,1:196], as.numeric)
trade2009 <- read.csv("trade2009.csv", header=T, row.names=1)
trade2009[,1:196] <- lapply(trade2009[,1:196], as.numeric)

dis2009 <- read.csv("dis2009.csv", header=T, row.names=1)
dis2009[,1:196] <- lapply(dis2009[,1:196], as.numeric)


# For node attributes (GDP, democracy, military expenditure) of 2009
load("vcov3.RData")
attri2009 <- subset(vcov3, year==2009,
                    select=state: log_military)

net2009 <- network(as.matrix(Edge2009), directed = T)

set.vertex.attribute(net2009,names(attri2009),attri2009)

set.network.attribute(net2009,"ally2009",as.matrix(ally2009))

set.network.attribute(net2009,"FDI2009",as.matrix(FDI2009))
set.network.attribute(net2009,"trade2009",as.matrix(trade2009))

set.network.attribute(net2009,"dis2009",as.matrix(dis2009))



##### Create list
networks <- list()# create network object
networks[[1]] <- net2010
networks[[2]] <- net2011
networks[[3]] <- net2012
networks[[4]] <- net2013
networks[[5]] <- net2014
networks[[6]] <- net2015
networks[[7]] <- net2016
networks[[8]] <- net2017
networks[[9]] <- net2018

# Set network attributes
network.vertex.names(x=net2009) <- attri2009$state
network.vertex.names(x=net2010) <- attri2010$state
network.vertex.names(x=net2011) <- attri2011$state
network.vertex.names(x=net2012) <- attri2012$state
network.vertex.names(x=net2013) <- attri2013$state
network.vertex.names(x=net2014) <- attri2014$state
network.vertex.names(x=net2015) <- attri2015$state
network.vertex.names(x=net2016) <- attri2016$state
network.vertex.names(x=net2017) <- attri2017$state
network.vertex.names(x=net2018) <- attri2018$state


LagNet<-list()
LagNet[[1]] <- net2009      # add network to the list
LagNet[[2]] <- net2010
LagNet[[3]] <- net2011
LagNet[[4]] <- net2012
LagNet[[5]] <- net2013
LagNet[[6]] <- net2014
LagNet[[7]] <- net2015
LagNet[[8]] <- net2016
LagNet[[9]] <- net2017


FDI2009<-network(as.matrix(FDI2009),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)

FDI2010<-network(as.matrix(FDI2010),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)


FDI2011<-network(as.matrix(FDI2011),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)

FDI2012<-network(as.matrix(FDI2012),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)

FDI2013<-network(as.matrix(FDI2013),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)

FDI2014<-network(as.matrix(FDI2014),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)

FDI2015<-network(as.matrix(FDI2015),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)

FDI2016<-network(as.matrix(FDI2016),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)

FDI2017<-network(as.matrix(FDI2017),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)

FDI2018<-network(as.matrix(FDI2018),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight", directed = T)


# Set network attributes
network.vertex.names(x=FDI2009) <- attri2009$state
network.vertex.names(x=FDI2010) <- attri2010$state
network.vertex.names(x=FDI2011) <- attri2011$state
network.vertex.names(x=FDI2012) <- attri2012$state
network.vertex.names(x=FDI2013) <- attri2013$state
network.vertex.names(x=FDI2014) <- attri2014$state
network.vertex.names(x=FDI2015) <- attri2015$state
network.vertex.names(x=FDI2016) <- attri2016$state
network.vertex.names(x=FDI2017) <- attri2017$state
network.vertex.names(x=FDI2018) <- attri2018$state


FDINet<-list()

FDINet[[1]] <- FDI2009
FDINet[[2]] <- FDI2010
FDINet[[3]] <- FDI2011
FDINet[[4]] <- FDI2012
FDINet[[5]] <- FDI2013
FDINet[[6]] <- FDI2014
FDINet[[7]] <- FDI2015
FDINet[[8]] <- FDI2016
FDINet[[9]] <- FDI2017

# For ally
ally2009<-network(as.matrix(ally2009), directed = F)
ally2010<-network(as.matrix(ally2010), directed = F)
ally2011<-network(as.matrix(ally2011), directed = F)
ally2012<-network(as.matrix(ally2012), directed = F)
ally2013<-network(as.matrix(ally2013), directed = F)
ally2014<-network(as.matrix(ally2014), directed = F)
ally2015<-network(as.matrix(ally2015), directed = F)
ally2016<-network(as.matrix(ally2016), directed = F)
ally2017<-network(as.matrix(ally2017), directed = F)
ally2018<-network(as.matrix(ally2018), directed = F)

# Set network attributes
network.vertex.names(x=ally2009) <- attri2009$state
network.vertex.names(x=ally2010) <- attri2010$state
network.vertex.names(x=ally2011) <- attri2011$state
network.vertex.names(x=ally2012) <- attri2012$state
network.vertex.names(x=ally2013) <- attri2013$state
network.vertex.names(x=ally2014) <- attri2014$state
network.vertex.names(x=ally2015) <- attri2015$state
network.vertex.names(x=ally2016) <- attri2016$state
network.vertex.names(x=ally2017) <- attri2017$state
network.vertex.names(x=ally2018) <- attri2018$state




allyNet<-list()
allyNet[[1]] <-ally2010
allyNet[[2]] <-ally2011
allyNet[[3]] <-ally2012
allyNet[[4]] <-ally2013
allyNet[[5]] <-ally2014
allyNet[[6]] <-ally2015
allyNet[[7]] <-ally2016
allyNet[[8]] <-ally2017
allyNet[[9]] <-ally2018


# For trade
trade2009<-network(as.matrix(trade2009),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)

trade2010<-network(as.matrix(trade2010),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)


trade2011<-network(as.matrix(trade2011),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)

trade2012<-network(as.matrix(trade2012),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)

trade2013<-network(as.matrix(trade2013),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)

trade2014<-network(as.matrix(trade2014),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)

trade2015<-network(as.matrix(trade2015),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)

trade2016<-network(as.matrix(trade2016),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)

trade2017<-network(as.matrix(trade2017),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)

trade2018<-network(as.matrix(trade2018),matrix.type="a", 
                   ignore.eval=FALSE,names.eval="weight_trade", directed = T)


# Set network attributes
network.vertex.names(x=trade2009) <- attri2009$state
network.vertex.names(x=trade2010) <- attri2010$state
network.vertex.names(x=trade2011) <- attri2011$state
network.vertex.names(x=trade2012) <- attri2012$state
network.vertex.names(x=trade2013) <- attri2013$state
network.vertex.names(x=trade2014) <- attri2014$state
network.vertex.names(x=trade2015) <- attri2015$state
network.vertex.names(x=trade2016) <- attri2016$state
network.vertex.names(x=trade2017) <- attri2017$state
network.vertex.names(x=trade2018) <- attri2018$state

# 
tradeNet<-list()

tradeNet[[1]] <- trade2009
tradeNet[[2]] <- trade2010
tradeNet[[3]] <- trade2011
tradeNet[[4]] <- trade2012
tradeNet[[5]] <- trade2013
tradeNet[[6]] <- trade2014
tradeNet[[7]] <- trade2015
tradeNet[[8]] <- trade2016
tradeNet[[9]] <- trade2017

# For distance
dis2009<-network(as.matrix(dis2009),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

dis2010<-network(as.matrix(dis2010),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)


dis2011<-network(as.matrix(dis2011),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

dis2012<-network(as.matrix(dis2012),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

dis2013<-network(as.matrix(dis2013),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

dis2014<-network(as.matrix(dis2014),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

dis2015<-network(as.matrix(dis2015),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

dis2016<-network(as.matrix(dis2016),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

dis2017<-network(as.matrix(dis2017),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

dis2018<-network(as.matrix(dis2018),matrix.type="a", 
                 ignore.eval=FALSE,names.eval="weight_distance", directed = T)

# Set network attributes
network.vertex.names(x=dis2009) <- attri2009$state
network.vertex.names(x=dis2010) <- attri2010$state
network.vertex.names(x=dis2011) <- attri2011$state
network.vertex.names(x=dis2012) <- attri2012$state
network.vertex.names(x=dis2013) <- attri2013$state
network.vertex.names(x=dis2014) <- attri2014$state
network.vertex.names(x=dis2015) <- attri2015$state
network.vertex.names(x=dis2016) <- attri2016$state
network.vertex.names(x=dis2017) <- attri2017$state
network.vertex.names(x=dis2018) <- attri2018$state

disNet<-list()
disNet[[1]] <-dis2010
disNet[[2]] <-dis2011
disNet[[3]] <-dis2012
disNet[[4]] <-dis2013
disNet[[5]] <-dis2014
disNet[[6]] <-dis2015
disNet[[7]] <-dis2016
disNet[[8]] <-dis2017
disNet[[9]] <-dis2018


# set overall network attributes

set.vertex.attribute(networks[[1]],"attri", attri2010)

set.vertex.attribute(networks[[2]],"attri", attri2011)

set.vertex.attribute(networks[[3]],"attri", attri2012)
set.vertex.attribute(networks[[4]],"attri", attri2013)
set.vertex.attribute(networks[[5]],"attri", attri2014)
set.vertex.attribute(networks[[6]],"attri", attri2015)
set.vertex.attribute(networks[[7]],"attri", attri2016)
set.vertex.attribute(networks[[8]],"attri", attri2017)
set.vertex.attribute(networks[[9]],"attri", attri2018)


# Test the model====

#  Model 1 in Table 2 ====
set.seed(1234)

fit3 <- btergm(networks ~
                 edgecov(FDINet, attrname="weight")
               +edgecov(LagNet)
               +edgecov(tradeNet,attrname="weight_trade")
               +edgecov(disNet, attrname="weight_distance")
      
               +edgecov(allyNet)
               
               + edges+ istar(2) + ostar(2) + mutual + ttriple + ctriple
               
               +nodeicov("log_military")
               +nodeicov("mean_Democracy")
               +nodeocov("log_GDP")
               
               +absdiff("mean_Democracy")
               +absdiff("log_GDP")
               +absdiff("log_military")
               +isolates, 
               R = 1000)
summary(fit3)
# Note:"Boot mean", not "Estimate", is used in Table 2 of the paper


# This time, I do not include network features (Model 3 in Table 2)====
set.seed(1234)
fit <- btergm(networks ~
                 edgecov(FDINet, attrname="weight")
               +edgecov(LagNet)
              +edgecov(tradeNet,attrname="weight_trade")
              +edgecov(disNet, attrname="weight_distance")
               
               +edgecov(allyNet)
               
               + edges
               
               +nodeicov("log_military")
               +nodeicov("mean_Democracy")
               +nodeocov("log_GDP")
               
               +absdiff("mean_Democracy")
               +absdiff("log_GDP")
               +absdiff("log_military")
               +isolates, 
               R = 1000)
summary(fit)
# Note:"Boot mean", not "Estimate", is used in Table 2 of the paper




################################
# Get the FDI graph for model 1 ( fit3 model): Upper left plot of Figure 3 ====

alldyads2 <- edgeprob(fit3)

# Fit the model manually to calculate confidence intervals

model <- glm(probability ~ edgecov.weight, family = "binomial", data = alldyads2)

# Create a data frame for predictions
new_data <- alldyads2 %>%
  mutate(pred = predict(model, newdata = ., type = "response"),
         se = predict(model, newdata = ., type = "link", se.fit = TRUE)$se.fit) %>%
  mutate(lower = pred - 1.96 * se,  # Lower bound of 95% CI
         upper = pred + 1.96 * se)  # Upper bound of 95% CI

# Plot with ggplot2
ggplot(new_data, aes(x = edgecov.weight, y = pred)) +
  geom_line(color = "blue") +                                # Main line
  geom_line(aes(y = lower), linetype = "dashed", color = "black") +  # Lower bound
  geom_line(aes(y = upper), linetype = "dashed", color = "black") +  # Upper bound
  labs(y = "Probability of Cosponsorship", x = "FDI (log)",caption = 
         "The effect of FDI of Model 1") +
  theme_minimal() +
  theme(
    axis.title.x = element_text(size = 24),  # X-axis label size
    axis.title.y = element_text(size = 24),  # Y-axis label size
    axis.text.x = element_text(size = 24),   # X-axis numbers (ticks) size
    axis.text.y = element_text(size = 24),    # Y-axis numbers (ticks) size
    plot.caption = element_text(size = 24, face = "bold", hjust = 0.5)  # Bigger, centered caption
  )



# Get the trade graph for model 1 ( fit3 model) Upper right plot of Figure 3 ====
model <- glm(probability ~ edgecov.weight_trade, family = "binomial", data = alldyads2)

# Create a data frame for predictions
new_data <- alldyads2 %>%
  mutate(pred = predict(model, newdata = ., type = "response"),
         se = predict(model, newdata = ., type = "link", se.fit = TRUE)$se.fit) %>%
  mutate(lower = pred - 1.96 * se,  # Lower bound of 95% CI
         upper = pred + 1.96 * se)  # Upper bound of 95% CI

# Plot with ggplot2
ggplot(new_data, aes(x = edgecov.weight_trade, y = pred)) +
  geom_line(color = "blue") +                                # Main line
  geom_line(aes(y = lower), linetype = "dashed", color = "black") +  # Lower bound
  geom_line(aes(y = upper), linetype = "dashed", color = "black") +  # Upper bound
  labs(y = "Probability of Cosponsorship", x = "Trade (log)",caption = 
         "The effect of trade of Model 1") +
  theme_minimal()+
  theme(
    axis.title.x = element_text(size = 24),  # X-axis label size
    axis.title.y = element_text(size = 24),  # Y-axis label size
    axis.text.x = element_text(size = 24),   # X-axis numbers (ticks) size
    axis.text.y = element_text(size = 24),    # Y-axis numbers (ticks) size
    plot.caption = element_text(size = 24, face = "bold", hjust = 0.5)  # Bigger, centered caption
  )




# To get the change statistics
load("FDI.RData")

df$log_FDI_a <- log(df$FDI_a+1)

df<-df[complete.cases(df), ] 

# get rid of some countries not need

df <- df[!df$country1_1=="AIA" & 
           !df$country1_1=="ABW" & 
           !df$country1_1=="BMU" & 
           !df$country1_1=="FLK" & 
           !df$country1_1=="GRL" & 
           !df$country1_1=="GLP" & 
           !df$country1_1=="GUF" & 
           !df$country1_1=="MTQ" & 
           !df$country1_1=="MSR" & 
           !df$country1_1=="MAF" & 
           !df$country1_1=="CUW" & 
           !df$country1_1=="SPM" & 
           !df$country1_1=="CYM" & 
           !df$country1_1=="REU" & 
           !df$country1_1=="COK" & 
           !df$country1_1=="FRO" & 
           !df$country1_1=="GIB" & 
           !df$country1_1=="GUM" & 
           !df$country1_1=="NCL" & 
           !df$country1_1=="SHN" & 
           !df$country1_1=="ASM" &
           !df$country1_1=="HKG" &
           !df$country1_1=="MAC" &
           !df$country1_1=="TWN" &
           !df$country1_1=="PYF",] 



df <- df[!df$country2_2=="AIA" & 
           !df$country2_2=="ABW" & 
           !df$country2_2=="BMU" & 
           !df$country2_2=="FLK" & 
           !df$country2_2=="GRL" & 
           !df$country2_2=="GLP" & 
           !df$country2_2=="GUF" & 
           !df$country2_2=="MTQ" & 
           !df$country2_2=="MSR" & 
           !df$country2_2=="MAF" & 
           !df$country2_2=="CUW" & 
           !df$country2_2=="SPM" & 
           !df$country2_2=="CYM" & 
           !df$country2_2=="REU" & 
           !df$country2_2=="COK" & 
           !df$country2_2=="FRO" & 
           !df$country2_2=="GIB" & 
           !df$country2_2=="GUM" & 
           !df$country2_2=="NCL" & 
           !df$country2_2=="SHN" & 
           !df$country2_2=="ASM" &
           !df$country2_2=="HKG" &
           !df$country2_2=="MAC" &
           !df$country2_2=="TWN" &
           !df$country2_2=="PYF",] 


# Explore how FDI affects the DV networks
FDI2018 <- df

change <- sd(FDI2018$log_FDI_a)
# 8.833


plogis(change*coef(fit3)[['edgecov.weight']])
# 10-29-2024
# [1] 0.5105196


# get the trade IV 
load("trade.RData")

df<-df3

df$log_trade_a <- log(df$trade+1)

df<-df[complete.cases(df), ] 

change <- sd(df$log_trade_a)
# 6.906


plogis(change*coef(fit3)[['edgecov.weight_trade']])
# 10-29-2024
# [1] 0.5333441



# Replication for Table 3 ====

# plogis(coef(fit3)[['edges']])
# [1] 0.02468658
# Edges are not shown in Table 3 to save spacex

# For istar: popularity
plogis(coef(fit3)[['istar2']])

# [1] 0.5109468


#sociability
plogis(coef(fit3)[['ostar2']])
# [1] 0.5083035

# mutuality
plogis(coef(fit3)[['mutual']])
# [1] 0.6013083

# transitivity
plogis(coef(fit3)[['ttriple']])
# [1] 0.5158833

# cyclicity
plogis(coef(fit3)[['ctriple']])
# [1] 0.4650287


# Replication for Table 4 ====

# # two sd of FDI
# change <- sd(FDI2018$log_FDI_a)*2
# # 17.666
# 


# For Network control effect
#Lag co-sponsor network
plogis(coef(fit3)[['edgecov.LagNet[[i]]']])

# [1] 0.8178594

# Ally network
plogis(coef(fit3)[['edgecov.allyNet[[i]]']])
# [1] 0.6967159

#Trade network
plogis(coef(fit3)[['edgecov.weight_trade']])
# [1] 0.5048353

#Distance
plogis(coef(fit3)[['edgecov.weight_distance']])
# [1] 0.4850566
# Note: negative sign is from the coefficient

# for homophily effects

# mean democracy

load("vcov3.RData")
attri2018 <- subset(vcov3,
                    select=state: log_military)

change <- sd(abs(apply(combn(attri2018$mean_Democracy,2), 2, diff)))
# [1] 0.2461481

# Democracy
plogis(change*coef(fit3)[['absdiff.mean_Democracy']])
# [1] 0.4465907
# Note: negative sign is from the coefficient; the same below for negative signs in tables

# log GDP
change <- sd(abs(apply(combn(attri2018$log_GDP,2), 2, diff)))
#5.0608

plogis(change*coef(fit3)[['absdiff.log_GDP']])
# [1] 0.4605549


# log military
change <- sd(abs(apply(combn(attri2018$log_military,2), 2, diff)))
# 8.9674

plogis(change*coef(fit3)[['absdiff.log_military']])
# [1] 0.5189836

# For node effects
# Log military
plogis(coef(fit3)[['nodeicov.log_military']]*(sd(attri2018$log_military)))
# [1] 0.5192986

# Democracy
plogis(coef(fit3)[['nodeicov.mean_Democracy']]*(sd(attri2018$mean_Democracy)))
# [1] 0.4683183

# log GDP
plogis(coef(fit3)[['nodeocov.log_GDP']]*(sd(attri2018$log_GDP)))
# [1] 0.4507698



### Interpret your model for specific states (10/29/24 updated)
# interpret(modelnamehere, i="n19", j="n20", type="tie") 
interpret(fit3, type = "tie", i = 32, j = 33, t =9)
# 0.01795196 

# Table 5 ====
interpret(fit3, type = "dyad", i = 32, j = 33, t = 9)
#          j->i = 0     j->i = 1
# i->j = 0 0.97345919 0.0085112021
# i->j = 1 0.01779495 0.0002346551

# Table 6 ====
interpret(fit3, type = "node", i = 32, j = c(33,85), t = 9)

#                probability Receiver 33 Receiver 85
# Sender 32 0.9596185520           0           0
# Sender 32 0.0004236051           1           1
# Sender 32 0.0175419459           1           0
# Sender 32 0.0224158970           0           1



# Online Appendix ====
# goodness of fit for Figure A1 ====

gof2<-gof(fit3, control=control.gof.ergm(nsim=1000))
plot(gof2)

# png("myplot2.png")




# goodness of fit for Figure A3 ====

gof2<-gof(fit, control=control.gof.ergm(nsim=1000))
plot(gof2)

# png("myplot2.png")


# End of script

































